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ABSTRACT 


This work traces the response of a granular material via the 
Ten Coefficient Truesdell rate-type constituitive model into the 
simplest meaningful loading: the triaxial test configuration. A 
functional relation has been posed for computing the rather peculiar 
relation between average applied stress and average porosity. Using 
that relation an attack has been mounted on the dilemna that 
exists between dynamic and constitutive use of the pressure variable; 
that is relating dynamic pressure, thermodynamic pressure, stress 
deviator and higher stress invariants. The resolution was as a 
linear superposition with a one-way feedback, in that while the 
dynamic component could not effect the constituitive component, 
the converse was not true since density appears in the momentum 
transport relation. 

There were two stages of preparation for this effort. The first 
was to use the PHOENICS Satellite to create the boundary conditions 
simulating the triaxial test. Two opposed parallel disks are joined 
with a thin cylindrical membrane with the granular material inside and 
a confining pressure outside. The membrane condition is modeled as 
a circular tension ring active on slices parallel to the disks; ie 
"slabs". Use was made of the PHOENICS provision for porosity greater 
than unity to represent buldging. The second preparation was for the 
granular constitution, which was coded into the Ground portion of the 
PHOENICS package. Provision was made by an ancillary computation to 
test for yield condition attempt to locate shear banding. Several non- 
granular cases were built to test the boundary conditions independent 
of the constitution and the low velocity response of PHOENICS and 
to disarm inherent simple fluid prejudice in the code. Some results 
were obtained and a covey of suggestions for parameter and computational 
control were generated. 
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INTRODUCTION 


This study has been directed at supporting an established and 
ongoing preparation for a triaxial test of spherical glass beads 
to be operated in the low-gravity environment of the Space Shuttle 
on orbit. The accronym MGM, mechanics of granular media, has 
been coined by the principal investigators. Dr. N.C. Costes and 
Prof. S. Sture, as the descriptor of this activity, and it is into 
that general arena that this effort has been cast. 

The first phases of this work were performed in the summer of 1984 
and reported in the image of this volume as Section XIV; "The 
Compatibility of an Existing CFD Code with a Broader Class of 
Constitutions" . The results and nomenclature of that report are 
to some extent assumed as background, especially the portions 
directed at the flexibilities built into the PHOENICS CFD Code. 
Reference will be made on occasion to groupings of variables 
unique to that report and to the general forms of transport density 
and pressure. 

To some extent the preliminary planning has been preoccupied with 
some of the aspects of getting the experiment into position for loading 
These are all fairly extensive problems in their own right. Some of 
them are: 1) consolidation by launch acceleration with vibration, 

2) membrane leakage during drawdown during confinement, 3) consequences 
of granule friction on the end platens during preload cycling, 4) the 
visibility of shear-bands with non-orientable granules ( spheres ) , 
and 5) the prejudice of a failure surface to a noncylindrical initial 
boundary. 

Granular mechanics with its more complicated constitutive behavoir 
has had to wait longer for analytic treatment than simple fluids and 
simple solids. The general realm of constitutive types is shown as 
Figure 1. In looking to associate a complicated behavoir with "fluid" 
or "solid", the back and forth coupling lines on the figure are very 
illuminating. The point being that a step up the tree toward 
generality from either simple fluid or simple solid branches is 
coupled to both. Hence it is perhaps not surprising that a flexible 
scheme for CFD might be ammenable to computational MGM. 

The governing parameter is the extent of intergranular stresses. In 
gravity currents of dust laden air or silt laden water 7 the granular 
material supplements the density but can be ignored in the constitution 
In a fluidized bed the granular matrix remains stationary or slowly 
circulating, but once levitated the intergranular stresses are low, 
likewise in the liquid counterpart percolation. Alternately in the 
very similar "pore-water" problem the active element is the granular 
matrix with the liquid adding a deviator hence the conditions; drained/ 
undrained. Recently there has been a flurry of activity using kinetic 
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theory to model high-speed granular flows, obviously a class cf 
flows in which the effects of the environmental fluid are ignored. 

A CFD package like PHOENICS includes features for multiphase flows 
to make modeling of granular-laden fluids possible, but of course 
allow no automatic mechanism for intergranular stress computation. 
Further Phoenics provides a constant named DARCY which builds a linear 
pressure gradient on volume flow rate. Herein quite a different 
approach is taken in which the constitution is coded into the exchange 
coefficients and functionally related to gradient of velocity -a^d 
the density as well as a "pressure” adjustment which adds intergranular 
stresses to the dynamic component. 

This study is then at the intersection of a contemporary CFD package 
taken to the outer limit of its flexibility and numerical MGM 
takenas its most comprehensive constitutive description. The emphasis 
is on matching feature for feature rather than creating a massive 
supplement to the PHOENICS code. 


/ 



FIGURE 1 Constitutive Realms 
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OBJECTIVES 


The primary objective of this study was to produce a cataloging 
of constitutive and boundary features that must be included in 
numerical models wether the destination code be PHOENICS or some 
alternative. While "PHOENICS thought" has been used to give some 
explicit definitions to the features, the features themselves are 
generic. Hence the object sought was the simplest model that retained 
the full features of the rationally derived MSD procedure while 
applying only to the triaxial configuration. Perhaps the greatest 
accomplishment in this regard is to arrive at specific conclusions 
for handling the boundary conditions in a way that is ammenable 
to the PHOENICS code, MICROFEM, or a custom code and more importantly 
to include the discoveries of preliminary testing. Further 
objectivity was sought by making explicit the offset terms in the 
Jaumann stress flux for this geometry and loading. 

Specifically the intention herein was to approach the goal of 
predicting specific loading responses from the rational direction, 
carefully highlighting the simplifications to permit backtracking should 
comparison of the model with experimental results ultimately show 
it to be inadequate due to over simplification. Since this study 
is preliminary to extensive testing it was not desireable to over 
simplify just to obtain "instant gratification" as analytic solutions. 

Another objective has been to emphasize the need for preparation of a 
numerical approach that explicitly generates fields of "signed" 
deviatoric stress power and makes them available to the evolving 
calculation. It is also desireable to determine the inhibition 
or enhancement parameter for shear bands based on granule size and 
shape, viz large spherical granules. 
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THE PRESSURE-DENSITY DILEMMA 


The goal is to relate pressure in the PHOENICS sense with the elements 
of tr ©* and to trace the flow of influence of density into the 
PHOENICS pressure both inherently and user coded. The goal is extended 
to include the realm of the triaxial test and to show what measures 
are necessary to circumvent PHOENICS provisions and to implement 
the full constitution of a granular media. 

Seveness to the ancient Hebrews indicated completion of a cycle and 
so it is with PHOENICS; six adjacent neighbors in space and one in time. 
Figure 2 shows PHOENICS pressure as the pressure component with a 
dynamic origin, P DYN , and in obedience to Pascal's Principle it has 

equal components and lies along the hydrostatic line. The source term 
for momentum transport equations ( Vj^ ) at dynamic equilibrium is 

the gradient of pressure; hence: 

- V/W * f w) * ^ 

The appearance of density,^, and transport density , are the 

coupling between the consCitutive equation and the computation of pressure. 
Figure 4 illustrates a "Ring Model" of computation to highlight the 
numerical symbiosis of pressure, momentum and mass. If thermodynamic 
pressure is defined in terms of the stress tensor as; 

P = 1/3 tr 

TH 

and further if thermodynamic pressure is considered to be a linear 
superposition of a constitutive and a dynamic part; 

P = P + P 

TH CON DYN 

which would still appear along the hydrostatic axis for some types 
of constitutions viz; ideal gas. The linear assumption is of course 
not correct in general since the components have separate functional 
dependencies on density. In fact with a granular media behavoir 
changes during loading must be included and so tests have to be made 
for yield or rupture to determine the local functional relation between 
pressure and density. As a further complication, granular materials 
require a rate-type constitutive model, ( all the references but #2 deal 
with this necessity ) . Otherwise the time dependency in the thermodynamic 
pressure could be isolated in the dynamic part, which would be very 
convenient in the CFD environment. The lure of incremental models is 
that time can be "eliminated" by multiplying the constitutive equation 
by Dt. But at best the incremental formulation introduces coordinate 
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system problems more difficult for CFD than the gradient of velocity 
components which are necessary for rate-type models with stress flux. 

**AAbtb* - 

c- = L ] "X * [ 1 1> ♦ [ 3 «r 

I 

cnw = tyv<zr 

t 

<cr - c 


er =■ c 


To see what reductions can be made with a triaxial configuration , 
a glance at Figure 2 shows the profile of the triaxial test in 
normal stress space ( Ox. = 0^ ) , the hydrostatic condition ( <X= <%= o%) 

which is a subset, and the stress deviator ( (T[+ o£+ C ) which 

has the same direction, but less restriction. Consider first the 
normal stress-normal strain relationship as shown in Figure 5; an 
accepted shape in soils testing. The very important "knee region" 
shows the difficult dependence on initial density; double-valued 
for initially dense conditions and for large axial strain approaching 
the critical state asymptotically; for which it is impossible to 
distinguish between initially loose and initially dense. It is 
analytically convenient to introduce without any rational justification 
( and apparently without precedence ) the three parameter relation: 

SR/ SR = tanh( ae ) + b( 1 - [tanh (ae - c) ] ) 

C 

The triaxial case gives: p = 1/3 tr Cf = 1/3 0^ ( SR + 2 ) 

with: SR = 0 £/ Oi . The slope of the lower portion 

of Figure 5 is a sort of 3-D Poisson's Ratio: isC-i^O 

it is notable that the zero of the ordinate represents : tr D = C. 

While the loose extreme compresses monotonically , the dense extreme 
exhibits an initial compression followed by expansion to critical 
state. This change of behavoir can prove to be very destabilizing 
to a numerical scheme. 

The classical description of the critical state is given as; 

P v — v 

c / P = exp( ° c ) . , . , .... 

' o in which specific 

K 

volume, the inverse of density, has been used and P = 0^ - 0^ 


A Uur&rr to £j dUJ 




J ^ *•••£/ CTTi&C '/> L 
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so that tr (T = P -1/3 0“^ ( 1 + 2*SR ) . . . .O^ - confining pressure. 

comprehensive rate-type model, this behavoir is nested in the first 
term of the denominator of the pressure expression: 

p = - ^ 

and if all the thermodynamic time dependence were in the dynamic 
component of pressure; 

P= - P M / ( • • • 

The exponent y represents the nonclassical coupling between pressure 
and density. One of the difficulties of using this pressure model 
is that the deviatoric stress power can change sign in some regions 
of the test and the resulting loading/unloading behavoir must be 
traced by a supplementary computation in a CFD environment. 

Isochoric plastic flow is one ostensibly at constant volume in which 

tr D = 0 hence = 0, a flow with identically constant density. 

For this case the full constitutive relation reduces to: 

O t 

- O'* ~ v/o* + O* 'V/ 

and by definition the deviatoric stress rate tensor is : 

c% - C PcX- OsOih .03ED +M 2 D 

The velocity vector has ifefr gradient distributed in the standard way 
into deformation rate and rotation tensors : 

VV = D + W 


the assistance of the latter has already been invoked in making an 
objective quantity from the substantial derivative of stress. The 
deformation rate is inherent in the rate-type constitution and of 
course appears there in the stress power with deviatoric and normal 
components and in the ytold condition also in the stress power. 

<o~ O ~ + P "ta. P> 

Macroscopically thermodynamics will place a condition on the integral 

/ Lit JCWOJv]j-t 

which represents the total energy expenditure on the test material. 

o = x># + X- 
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FIGURE 2 


NORMAL STRESS SPACE 


Failure comes in the forms of rupture and yield. Figure 2 shows 
on a triaxial plane the lack of symmetry of the rupture boundaries 
about the hydrostatic axis ( "compression/ extension failure" ) * 

The other surface forming the failure boundary is given by the yield 
cap. Following the procedure of MSD ( 78 - 84 ) the yield surface 
is established by a set of Q vectors which render the array G, singular 

D 2 

Q = [G ( Q f jo >] iP . . . Q = [ P 2 ] 

q 

p - stress deviator. . . 1/3 tr 0 

q - square root of the second invariant \/ 5 * 6 * q = M o c 

p tr D C 

IP - power tensor - [ ] 

tr <7* D* 


The density variation coming from the critical pressure-specific 
volume relation given above. An overly simple verbal description 
of this equation might run; the stress invariants are related to 
exponentially decaying transients scaled in the components of stress 
power. Aside from that description an algebraic result can be 
obtained if the reduced constitution ( isochoric) above is used; 

q 2 = ( 2 p c p + p 2 ) 

Mq and p^ are the classical parameters of critical state. This 
relation plots as an ellipse in q-p space and hence the failure 
envelop looks like a somewhat distorted icecream cone; as it appears 
in sections of the q-p- jO space at constant density. 









Digressing to the specifics of a triaxial configuration with 
full 0 symmetry; the stress deviator is: 

l/3tr Oi +20^ 

and hence the deviatoric stress is: 

[ 2/3(0 rj-Oj ] 




and the deformation rate is : 

d = ‘4 I 
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and the spin rate is : 


All these components then can be used with the exchange coefficients 
for the PHOENICS momentum transport equations; for the ith : 


o 

c* 


I " ’»* - 

f <k = o)4(v/c~* = cr*'W') 
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FIGURE 5 A RING MODEL FOR PHOENICS 




PRESSURE-DENSITY STATES 
FIGURE 3 


GRANULAR STRESS-STRAIN 
& POISSON 
FIGURE 4 





FIGURE 9 SLABS OF SPHERES 


PLATEN FRICTION 

(pseudo-bulging) 

LOCAL CURVATURE 


SHEAR BAND 
LOCATION/ DETECTION 


CONSOLIDATION 

HISTORY 


CONFINEMENT 

VS 

MEMBRANE ELASTICITY 



FIGURE 10 SPECIAL CONSIDERATIONS 
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FIGURE 7 

List' of PHOENICS FILES (KENFILES.TXT) 


# XXX 


Implementations 


2D POLAR PORUS PLUG 10X20X1 
UNIFORM INFLOW INTO ANNULUS 
FREE REGION AROUND PLUG 
3D CARTESIAN DIAMOND FIELD 
20X20X20 OBSTRUCTIONS 
UNIFORM INFLOW 

DEL SQR T = DT/Dt. .CONDUCTION 
TRANSIENT 20X20X1 2D NOFLOW 
DEL SQR Hi = KCl & DEL SQR H2 = KC2 
BUILD FILE.. LOW V IN KWF 
TRI AXIAL LOADING 

TRIAXIAL GEOMETRY WITH IDEAL GAS 
ENCHROACHING END CAP 
DRAINING UNDER VACUUM 
MASS REDISTRIBUTED OVER MULTIRUN 
MEMBRANE CONDITION/ LATERAL 
TRIAXIAL LOADING: HYPERBOLIC TANH 
PREYIELD: "EXCO" : GRANULAR 
(TRI AND TRY AGAIN ! ) YIELD 
DET G SINGULAR PORTION OF FIELD 
BENARD: LOW VELOCITY TEST WITH ONSET 

OF CONVECTION .-PLANE LAYER 10X10X10 
BENFLD INITIAL FIELD FOR INSTABILITY 
USER NAME KWF 

USER NAME KWF . . RUNS TALSAT.FTN, .EAR,. LOG 



FIGURE 8 POROSITY VARYING WITH SCALED PRESSURE 


1 KEN 

2 FRE 


3 

KWF 

4 

NCH 

5 

NET 

6 

TAL 

7 

TAG 

8 

TAGl 

9 

TAG 2 

10 

TAG 3 

11 

TAG 4 

12 

TRI 

13 

TRY 

14 

BEN 


15 KEN (.JOB) 

16 TAL (.JOB) 
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LATERAL & LONGITUDINAL BOUNDARIES 


Aside from dealing with a very difficult constitution internal to the 
cylindrical test region , the boundary conditions present formidable 
difficulties to a CFD package; Figures 2,9&10. The end platens are 
rigid and nominally frictionless; however preliminary testing has 
shown that if the granular material is glass spheres ( a necessity 
for a SS based experiment ), the friction is very tenuous. The 
consequences of friction to the boundary conditions and to numerical* 
computation are that an additional confinement is present on the 
upper and lower regions ( 1/3 platen diameter ? ) of the cylinder of 
test material. Typically this over-confinement causes premature, 
asymmetrical bulging. Further there will certainly be a disruption 
to the proper growth of the failure surfaces; viz. shear bands. 

The steadily advancing ( constant c Figure 4 ) upper platen must be 
treated with caution numerically, since a step displacement might cause 
a jump in local density and hence an irreversible catastrophic change 
in the local constitution. That is the parallel running computations 
for testing realms of behavoir might be triggered, before the density 
wave could diffuse. ... "wave and diffuse referring to the numerical and 
not physical uses of the terms". 



This advancing platen problem has been tackled in the same way as the 
confining membrane condition. In both a use of the porosity in the 
transport density, ^R' was implemented. The platen is replaced by 
an advancing wave of increasing porosity; ie the "fuzzy platen". For 

the bulging of the membrane it will be treated 
by allowing the circumferencial cells to have 
porosity in excess of one. Figure 8 shows a 
i simple model of the change of porosity with 

which is a pressure scaled to membrane stress, 
P/2MK. Since friction over the spherical 
surfaces changes the membrane tension exponentially 
( T/Tq = exp(-fe)), the lateral membrane effects 
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are small and a narrow slab ( say 10 D 0 ) can be isolated for analysis. 
This will give an artificial "aliasing" as the resolution jumps from 
D q to 10 *D q and since shear bands will be suppressed to vertical 
orientation on each slab. 


RESULTS & CONCLUSIONS 


The preliminary specifics for the numerical treatment of triaxial 
loading of a granular material in the PHOENICS environment have been 
laid to rest; constitution, boundary conditions and loading simulation; 
with the details given herein for slabwise, 0 independent modeling. 
Figure 6 shows some of the forays into the numerical realm to test 
various schemes for handling this configuration and loading. Special 
interest was placed on the treatment of dynamic and thermodynamic 
pressure and the correct decomposition of stress tensor into that 
pressure and the momentum exchange coefficient. The 0 dependent terms 
in Grad V, 0, etc, if the shear bands are to be predicted. That 
inclusion will be necessary but advisable only after the much simpler 
version is compared with test results. Those 0 dependent terms make 
a considerable change in the accounting and in the "objectivity offset", 
but the changes are straightforward in the current framework. 

A continuing doubt of the validity of PHOENICS at very low 
velocities and especially at small displacements, ( V dt ) persists. 
However the preparations and the lessons in coding will be useful in 
general and to any numerical approach to the triaxial loading of a 
granular material. A considerable portion of the summer's effort was 
spent pursuing minor avenues ( draining, percolation, preloading, etc ) 
that are not included in this report but which provide support to the 
project as a whole and all delt with various aspects of MGM. 
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